__author__ = 'Benjamin Weaver <benjamin.weaver@noirlab.edu>'
__version__ = '202100302'
__datasets__ = ['sdss_dr16']
__keywords__ = ['HowTo','spectra','query','visualization','SDSS']
Obtain SDSS spectra using the Astro Data Lab spectrum service, convert data to specutils objects as needed, and use Prospect to display the data.
Prospect is an interactive spectrum visualization service that is being actively used by the DESI project for visual inspection of commissioning and survey validation spectra right now. Rather than just a simple flux versus wavelength plot, the interactive spectrum display provides pan, zoom, known spectral lines, targeting and other catalog information, metadata, etc. The Bokeh library handles the low-level interactivity.
We can leverage the generic spectrum container objects provided by specutils to allow Prospect to plot data from other surveys. In this example, we are using SDSS DR16.
In the future we will extend this functionality to additional instruments. For example, we are excited by the possibility of using Prospect and the Astro Data Lab spectrum service for Gemini spectra. And, of course, once DESI data become public, both the spectrum service and this tool will be ready.
If you use this notebook for your published science, please acknowledge the following:
This cell performs all necessary imports for this notebook.
import os
import sys
import time
from getpass import getpass
#
# Numpy, astropy, etc.
#
import numpy as np
import astropy.units as u
from astropy.table import Table
from astropy.nddata import InverseVariance
from specutils import Spectrum1D, SpectrumCollection, SpectrumList
#
# Prospect imports
#
from prospect.specutils import read_spPlate, read_spZbest
from prospect.viewer import plotspectra
#
# Astro Data Lab imports
#
from dl import specClient as spec # primary spectral data client interface
from dl import storeClient as sc # needed to use virtual storage
from dl import queryClient as qc # needed to query Astro Data Lab catalogs
from dl import authClient as ac # needed for login authentication
start_time = time.time() # save start time of notebook run
Much of the functionality of spectrum services can be accessed without explicitly logging into Astro Data Lab (the service then uses an anonymous login). But some capabilities, for instance saving the results of your queries to your virtual storage space, require a login (i.e. you will need a registered user account).
If you need to log in to Astro Data Lab, uncomment the ac.login() command and respond according to the prompts. If you have previously logged into Astro Data Lab, this cell will simply print your active user name.
# ac.login(input("Enter user name: (+ENTER) "),getpass("Enter password: (+ENTER) "))
ac.whoAmI()
'demo01'
SDSS spectra are stored per-plate in spPlate files. These contain 640 spectra for the original SDSS spectrograph or 1000 spectra for the BOSS/eBOSS spectrograph. All spectra in an individual spPlate file have a common wavelength solution, so a single spPlate file can be represented by a Spectrum1D object. The Prospect visualization package can then take this Spectrum1D object as input data.
In the cells below, we will find spectra corresponding to SDSS plate 2955, but using the spectrum service instead of reading the spectra from a file. This starts by finding the coordinates of the plate center.
plate = 2955
q = f"SELECT ra, dec FROM sdss_dr16.platex WHERE plate = {plate}"
r = qc.query(sql=q, fmt='array')
r
array([235.4204 , 1.644156])
SDSS plates are circular, so we'll do a cone search around the plate center found above. This gets us the spectrum identifiers we need to retrieve the spectra themselves.
sdss_ids = spec.query(r[0], r[1], 1.5, fmt='array', out='', constraint=f'plate = {plate} ORDER BY specobjid LIMIT 50')
sdss_ids
array([3327035125891360768, 3327036500280895488, 3327037050036709376,
3327038424426244096, 3327038974182057984, 3327040073693685760,
3327040348571592704, 3327040623449499648, 3327040898327406592,
3327041173205313536, 3327041722961127424, 3327041997839034368,
3327042822472755200, 3327043097350662144, 3327043372228569088,
3327043921984382976, 3327044196862289920, 3327045021496010752,
3327045296373917696, 3327046121007638528, 3327046670763452416,
3327047495397173248, 3327048045152987136, 3327048320030894080,
3327048594908801024, 3327048869786707968, 3327049144664614912,
3327049419542521856, 3327049694420428800, 3327050244176242688,
3327051068809963520, 3327051343687870464, 3327052443199498240,
3327052992955312128, 3327053267833219072, 3327053817589032960,
3327054092466939904, 3327054642222753792, 3327054917100660736,
3327055191978567680, 3327055466856474624, 3327056016612288512,
3327056566368102400, 3327056841246009344, 3327057116123916288,
3327057391001823232, 3327057665879730176, 3327057940757637120,
3327058215635544064, 3327058490513451008])
Now we fetch the spectra. The option align=True directes the service to return all spectra as a single Spectrum1D object with a common wavelength axis for all spectra. Since all spectra on a single SDSS plate have a common wavelength axis anyway, this is a simple operation.
%%time
sdss = spec.getSpec(sdss_ids, fmt='Spectrum1D', align=True)
sdss
CPU times: user 27.5 ms, sys: 22.2 ms, total: 49.7 ms Wall time: 225 ms
<Spectrum1D(flux=<Quantity [[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]] erg / (Angstrom cm2 s)>, spectral_axis=<SpectralAxis [3802.77 , 3803.6448, 3804.522 , ..., 9217.223 , 9219.344 , 9221.464 ] Angstrom>, uncertainty=InverseVariance([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]))>
The Prospect viewer requires certain spectral metadata that are normally provided with SDSS spectra--i.e., these metadata are stored alongside the flux versus wavelength data. There are two categories of metadata:
In the cells below, we will extract the necessary metadata from the specobj and photoplate tables and reformat them slightly to match the appearance of these data as they would be extracted from SDSS files.
This query extracts the redshift and classification information from the specobj table.
q = "SELECT specobjid, class, subclass, z, zerr, rchi2diff, zwarning FROM sdss_dr16.specobj WHERE specobjid IN ({0}) ORDER BY specobjid".format(','.join(map(str, sdss_ids.tolist())))
sdss_z = qc.query(sql=q, fmt='table')
for c in ('specobjid', 'class', 'subclass', 'z', 'zerr', 'rchi2diff', 'zwarning'):
if c == 'zerr':
sdss_z.rename_column(c, 'Z_ERR')
else:
sdss_z.rename_column(c, c.upper())
sdss_z
| SPECOBJID | CLASS | SUBCLASS | Z | Z_ERR | RCHI2DIFF | ZWARNING |
|---|---|---|---|---|---|---|
| int64 | str6 | str11 | float64 | float64 | float64 | int64 |
| 3327035125891360768 | GALAXY | -- | 0.26219922 | 6.426468e-05 | 0.16153336 | 0 |
| 3327036500280895488 | GALAXY | -- | 0.21097003 | 5.7919853e-05 | 0.13921309 | 0 |
| 3327037050036709376 | STAR | K5 | -2.3494376e-06 | 1.4350197e-05 | 0.49699175 | 0 |
| 3327038424426244096 | GALAXY | -- | 0.2605426 | 6.229882e-05 | 0.22558784 | 0 |
| 3327038974182057984 | GALAXY | -- | 0.11276812 | 1.5994665e-05 | 0.510905 | 0 |
| 3327040073693685760 | GALAXY | -- | 0.06364931 | 2.8312488e-05 | 0.42436635 | 0 |
| 3327040348571592704 | GALAXY | -- | 0.20836118 | 4.189276e-05 | 1.1794078 | 0 |
| 3327040623449499648 | GALAXY | -- | 0.0 | 0.0 | 0.0 | 134 |
| 3327040898327406592 | GALAXY | -- | 0.20784536 | 4.1236628e-05 | 0.61972916 | 0 |
| 3327041173205313536 | QSO | -- | 0.04825218 | 0.00014616082 | 0.00018835068 | 4 |
| ... | ... | ... | ... | ... | ... | ... |
| 3327055466856474624 | GALAXY | -- | 0.16715993 | 3.366827e-05 | 0.5121331 | 0 |
| 3327056016612288512 | GALAXY | STARFORMING | 0.07140229 | 1.3711498e-05 | 0.6119038 | 0 |
| 3327056566368102400 | GALAXY | -- | 0.061552946 | 1.5951675e-05 | 1.0402132 | 0 |
| 3327056841246009344 | QSO | -- | 1.1373804 | 0.0006364495 | 0.06628287 | 0 |
| 3327057116123916288 | STAR | F9 | -0.00017936311 | 1.3192781e-05 | 0.8338379 | 0 |
| 3327057391001823232 | STAR | A0 | -0.00055759837 | 3.88684e-05 | 0.19958907 | 0 |
| 3327057665879730176 | STAR | M3 | 9.9042336e-05 | 7.070536e-05 | 0.092519104 | 0 |
| 3327057940757637120 | GALAXY | -- | 0.41268384 | 1.4011068 | 0.0028770566 | 5 |
| 3327058215635544064 | GALAXY | STARFORMING | 0.060996726 | 8.68608e-06 | 1.6069216 | 0 |
| 3327058490513451008 | QSO | BROADLINE | 0.82612693 | 0.00022074989 | 0.276842 | 0 |
This query extracts photometric targeting information from the photoplate table. The data obtained needs to be reformatted to match the "plugmap" metadata that is stored with SDSS spectrum files, including spPlate files.
q = """SELECT s.specobjid, s.targetobjid, s.ra, s.dec,
s.primtarget, s.sectarget,
s.boss_target1,
s.ancillary_target1, s.ancillary_target2,
s.eboss_target0, s.eboss_target1, s.eboss_target2,
p.u, p.g, p.r, p.i, p.z
FROM sdss_dr16.specobj AS s
LEFT JOIN sdss_dr16.photoplate AS p
ON s.bestobjid = p.objid
WHERE specobjid IN ({0})
ORDER BY specobjid""".format(','.join(map(str, sdss_ids.tolist())))
sdss_plugmap = qc.query(sql=q, fmt='table')
for c in ('ra', 'dec', 'primtarget', 'sectarget',
'boss_target1',
'ancillary_target1', 'ancillary_target2',
'eboss_target0', 'eboss_target1', 'eboss_target2'):
sdss_plugmap.rename_column(c, c.upper())
sdss_plugmap['OBJID'] = np.vstack((np.bitwise_and(sdss_plugmap['targetobjid'] >> 32, 2**16 - 1).data,
np.bitwise_and(sdss_plugmap['targetobjid'] >> 48, 2**11 - 1).data,
np.bitwise_and(sdss_plugmap['targetobjid'] >> 29, 2**3 - 1).data,
np.bitwise_and(sdss_plugmap['targetobjid'] >> 16, 2**12 - 1).data,
np.bitwise_and(sdss_plugmap['targetobjid'], 2**16 - 1).data)).T
sdss_plugmap['MAG'] = np.vstack((np.where(np.isnan(sdss_plugmap['u'].data), 0, sdss_plugmap['u'].data),
np.where(np.isnan(sdss_plugmap['g'].data), 0, sdss_plugmap['g'].data),
np.where(np.isnan(sdss_plugmap['r'].data), 0, sdss_plugmap['r'].data),
np.where(np.isnan(sdss_plugmap['i'].data), 0, sdss_plugmap['i'].data),
np.where(np.isnan(sdss_plugmap['z'].data), 0, sdss_plugmap['z'].data))).T
sdss.meta['plugmap'] = sdss_plugmap
sdss_plugmap
| specobjid | targetobjid | RA | DEC | PRIMTARGET | SECTARGET | BOSS_TARGET1 | ANCILLARY_TARGET1 | ANCILLARY_TARGET2 | EBOSS_TARGET0 | EBOSS_TARGET1 | EBOSS_TARGET2 | u | g | r | i | z | OBJID [5] | MAG [5] |
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| int64 | int64 | float64 | float64 | int64 | int64 | int64 | int64 | int64 | int64 | int64 | int64 | float64 | float64 | float64 | float64 | float64 | int64 | float64 |
| 3327035125891360768 | 11268994537292320 | 236.58816 | 0.91476191 | 96 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 21.303427 | 19.10137 | 17.512045 | 16.937172 | 16.563414 | 2327 .. 544 | 21.303427 .. 16.563414 |
| 3327036500280895488 | 11268994537161253 | 236.31498 | 0.83142165 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.57446 | 19.067663 | 17.930079 | 17.399467 | 16.985687 | 2327 .. 549 | 20.57446 .. 16.985687 |
| 3327037050036709376 | 11268994537095484 | 236.17657 | 0.91950364 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.861134 | 18.504723 | 17.763912 | 17.445564 | 17.245972 | 2327 .. 316 | 19.861134 .. 17.245972 |
| 3327038424426244096 | 11268994537292314 | 236.57992 | 0.84114607 | 32 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 21.99351 | 19.606033 | 18.029598 | 17.456076 | 17.053299 | 2327 .. 538 | 21.99351 .. 17.053299 |
| 3327038974182057984 | 11268994537292219 | 236.63067 | 0.85355411 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.648071 | 18.084349 | 17.2542 | 16.836885 | 16.52554 | 2327 .. 443 | 19.648071 .. 16.52554 |
| 3327040073693685760 | 11265262790246829 | 236.57247 | 1.3802039 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.00688 | 18.037868 | 17.137938 | 16.725918 | 16.394184 | 1458 .. 429 | 20.00688 .. 16.394184 |
| 3327040348571592704 | 11265262790312240 | 236.66305 | 1.4061234 | 96 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.11582 | 17.984798 | 16.497595 | 15.93136 | 15.503163 | 1458 .. 304 | 20.11582 .. 15.503163 |
| 3327040623449499648 | 11265262252261744 | 234.00554 | 1.1758617 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | -- | -- | -- | -- | -- | 1458 .. 368 | 0.0 .. 0.0 |
| 3327040898327406592 | 11268995074162967 | 236.57875 | 1.3024246 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 21.231564 | 18.893684 | 17.517618 | 16.979557 | 16.593437 | 2327 .. 279 | 21.231564 .. 16.593437 |
| 3327041173205313536 | 11265262790246975 | 236.60779 | 1.4217924 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 24.281784 | 21.341381 | 20.427032 | 20.001045 | 20.201952 | 1458 .. 575 | 24.281784 .. 20.201952 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 3327055466856474624 | 11265262790181186 | 236.35198 | 1.4705854 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.896185 | 18.988594 | 17.842579 | 17.310987 | 16.932745 | 1458 .. 322 | 20.896185 .. 16.932745 |
| 3327056016612288512 | 11265262790181282 | 236.43226 | 1.3820505 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.015394 | 17.45155 | 16.701193 | 16.29317 | 15.999505 | 1458 .. 418 | 19.015394 .. 15.999505 |
| 3327056566368102400 | 11268994536767942 | 235.44215 | 0.87831515 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.036299 | 18.109346 | 17.067467 | 16.503475 | 16.056887 | 2327 .. 454 | 20.036299 .. 16.056887 |
| 3327056841246009344 | 11268994536702418 | 235.30366 | 0.84224156 | 1048578 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 20.499498 | 20.08053 | 19.363874 | 18.998253 | 18.802431 | 2327 .. 466 | 20.499498 .. 18.802431 |
| 3327057116123916288 | 284707516645555 | 235.31638 | 0.59804179 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 19.503473 | 18.193771 | 17.647467 | 17.41578 | 17.298328 | 752 .. 179 | 19.503473 .. 17.298328 |
| 3327057391001823232 | 284707516645557 | 235.31813 | 0.43006452 | 0 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 19.230766 | 18.273815 | 17.969347 | 17.86463 | 17.81168 | 752 .. 181 | 19.230766 .. 17.81168 |
| 3327057665879730176 | 11268994536768101 | 235.41088 | 0.90478873 | 536877056 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 23.794682 | 22.39897 | 21.048717 | 20.003584 | 19.45056 | 2327 .. 613 | 23.794682 .. 19.45056 |
| 3327057940757637120 | 11268994536703488 | 235.24978 | 0.85021657 | 0 | 16 | 0 | 0 | 0 | 0 | 0 | 0 | -- | -- | -- | -- | -- | 2327 .. 1536 | 0.0 .. 0.0 |
| 3327058215635544064 | 11268994536767935 | 235.43123 | 0.86691538 | 64 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.057314 | 17.920666 | 17.439993 | 17.0802 | 16.939682 | 2327 .. 447 | 19.057314 .. 16.939682 |
| 3327058490513451008 | 11268994536767618 | 235.36244 | 0.8715034 | 1048578 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | 19.41755 | 19.03087 | 18.888653 | 18.939404 | 18.742466 | 2327 .. 130 | 19.41755 .. 18.742466 |
Once all the spectral data are assembled, the visualization is launched with a single plotspectra() function in the cell below. The primary input is the Spectrum1D object created above. Prospect expects a certain scaling on the input flux, so sdss.new_flux_unit() handles that (this changes the numerical scale factor on the flux unit, not the flux unit itself). Other inputs:
zcatalog: the redshift table created above.model: the model spectrum used by the SDSS pipeline to obtain the redshift and classification.notebook: instructs Prospect to use notebook display mode (as opposed to standalone HTML).title, model_from_zcat, etc.: these options control certain display options that are not directly related to the input data or that need to be switched off for SDSS data. A full description is available.Once the cell is run, it will take several seconds for the data to load and render. Then you will be able to pan and zoom on the primary flux versus wavelength display, page through different spectra, show known emission or absorption lines, etc. More details are below the display.
plotspectra(sdss.new_flux_unit(u.Unit('1e-17 erg/(cm**2 s Angstrom)')),
zcatalog=sdss_z,
model=(sdss.spectral_axis.value, sdss.meta['model'].to(u.Unit('1e-17 erg/(cm**2 s Angstrom)'))),
notebook=True, title=os.path.basename('SDSS Plate %d' % plate), model_from_zcat=False, with_coaddcam=False, mask_type='PRIMTARGET', with_thumb_tab=False, with_vi_widgets=False)
TARGETID: For SDSS data, this tuple of numbers corresponds to the photometric imaging run, photometric processing version, camera column ("camcol"), imaging field (within the run), and object ID within the field.mag_u, ...: These are the photometric magnitudes in SDSS u, g, r, i, z filters.SPECTYPE: The type of object, typically 'GALAXY', 'QSO', 'STAR'.SUBTYPE: A more specific classification of SPECTYPE, for example a 'STAR' may be 'K5' or a 'GALAXY' may be 'STARFORMING'.Z: The redshift.ZERR: The uncertainty on Z.ZWARN: If the processing pipeline produced errors, this value will be non-zero.For reference, and to compare timing, below is how one would fetch spectra and redshift metadata from file-based storage. First, obtain the path to the files in the Astro Data Lab file service.
spectro_redux = 'sdss_dr16://' + os.path.join('sdss', 'spectro', 'redux')
run2d = '26'
plate = 2955
mjd = '54562'
sdss_spectra = os.path.join(spectro_redux, run2d, str(plate), f'spPlate-{plate:d}-{mjd}.fits')
sdss_redshifts = os.path.join(spectro_redux, run2d, str(plate), f'spZbest-{plate:d}-{mjd}.fits')
print(sdss_spectra)
print(sdss_redshifts)
sdss_dr16://sdss/spectro/redux/26/2955/spPlate-2955-54562.fits sdss_dr16://sdss/spectro/redux/26/2955/spZbest-2955-54562.fits
Next, read the spPlate file. This automatically loads the targeting metadata.
%%time
sdss_file = read_spPlate(sc.get(sdss_spectra, mode='fileobj'), limit=50)
sdss_file
CPU times: user 349 ms, sys: 260 ms, total: 609 ms Wall time: 2.22 s
<Spectrum1D(flux=<Quantity [[11.0880165 , 11.084556 , 11.081101 , ..., 11.787075 ,
11.787052 , 11.787028 ],
[ 4.371129 , 4.36966 , 4.3681927 , ..., 2.2523484 ,
2.2523346 , 2.2523208 ],
[ 0.85558206, 0.8553027 , 0.8550237 , ..., 8.030891 ,
8.030855 , 8.03082 ],
...,
[12.411868 , 12.408562 , 12.405259 , ..., 25.452364 ,
25.452513 , 25.45266 ],
[-1.97795 , -1.9774021 , -1.976855 , ..., 15.674774 ,
15.674834 , 15.674895 ],
[ 0.6398605 , 0.63967246, 0.63948476, ..., -0.11592854,
-0.11592867, -0.1159288 ]] 1e-17 erg / (Angstrom cm2 s)>, spectral_axis=<SpectralAxis [3802.76948244, 3803.64520329, 3804.5211258 , ..., 9217.22098659,
9219.34357452, 9221.46665124] Angstrom>, uncertainty=InverseVariance([[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
...,
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.],
[0., 0., 0., ..., 0., 0., 0.]]))>
Read the spZbest file, which contains the redshift catalog.
sdss_z_file, sdss_model_file = read_spZbest(sc.get(sdss_redshifts, mode='fileobj'), limit=50)
Launch the visualization. Reading from files produces slightly different scaling on the flux values, but otherwise the options are the same as the example above.
plotspectra(sdss_file,
zcatalog=sdss_z_file,
model=(sdss_model_file.spectral_axis.value, sdss_model_file.flux.value,),
notebook=True, title=os.path.basename('SDSS Plate %d' % plate), model_from_zcat=False, with_coaddcam=False, mask_type='PRIMTARGET', with_thumb_tab=False, with_vi_widgets=False)
end_time = time.time()
print('Notebook run time: %.3g sec' % (end_time - start_time))
Notebook run time: 37.2 sec